library(dplyr)
library(magrittr)
library(tidyr)
library(ggplot2)

dat <- rio::import("data/Demo_book_country_3.dta")

dem <- dat %>%
  group_by(country_id) %>%
  arrange(year) %>%
  mutate(lp_lat_abst = zoo::na.locf(lp_lat_abst, na.rm = F),
         portdist_natural_km = zoo::na.locf(portdist_natural_km, na.rm = F)) %>%
  data.frame


dem %<>%
  filter(year >= 1789 & year <= 2019, europe_harbors_def == 0) %>%
  select(country_id, year, v2x_polyarchy_imp_100, 
         eur_pct_est_smooth, portdist_natural_km, lp_lat_abst)


regs_df <- data.frame(start_year = 1789:1990,
                      end_year = 1818:2019,
                      coef = NA,
                      se = NA,
                      t_value = NA,
                      N = NA,
                      r2 = NA, 
                      series = "lexical")

for (i in 1:nrow(regs_df)){
  subset_data <- dem[which(dem$year %in% c(regs_df$start_year[i]:regs_df$end_year[i])), ]
  reg <- miceadds::lm.cluster(v2x_polyarchy_imp_100 ~ eur_pct_est_smooth +
                                portdist_natural_km + lp_lat_abst + as.factor(year),
                              data = subset_data, cluster = "country_id")
  reg_sum <- summary(reg)
  regs_df$coef[i] <- reg_sum[2, 1]
  regs_df$se[i] <- reg_sum[2, 2]
  regs_df$t_value[i] <- reg_sum[2, 3]
  regs_df$N[i] <- nrow(reg$lm_res$model)
  regs_df$r2[i] <- summary(reg$lm_res)$r.squared
}

regs_df$year <- regs_df$end_year - 15

g <- regs_df %>%
  filter(series == "lexical") %>%
  ggplot(aes(x = year)) +
  geom_ribbon(aes(ymin = coef - se,
                  ymax = coef + se),
              color = 'lightgray',
              alpha = 0.25) +
  geom_point(aes(y = coef)) +
  scale_x_continuous(breaks = seq(1800, 2025, by = 25)) +
  labs(x = "Year",
       y = "Coefficient estimates for European Ancestry") +
  theme(panel.grid.minor = element_blank(),
        panel.background = element_blank(),
        panel.grid.major = element_blank(),
        panel.border = element_blank(),
        axis.line = element_line(color = "black"),
        plot.margin = margin(1, 1, 1, 0.25, "cm"))

## ggsave(plot = g, filename = "figure_11_2.png",
##        width = 11, height = 8.5)

ggsave(plot = g, filename = "output/figure_11_2.tiff",
       width = 11, height = 8.5, dpi = 300)
